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We compare weak and strong coupling theory of counterion-mediated electrostatic interactions 
between two asymmetrically charged plates with extensive Monte-Carlo simulations. Analytical 
results in both weak and strong coupling limits compare excellently with simulations in their re- 
spective regimes of validity. The system shows a surprisingly rich structure in terms of interactions 
between the surfaces as well as fundamental qualitative differences in behavior in the weak and the 
strong coupling limits. 



I. INTRODUCTION 

Stability and interactions in biological and soft-matter 
systems often depends on the underlying properties of 
electrostatic interactions Charges on macromolec- 
ular surfaces in aqueous environments, as in the case 
of membranes, self-assembled micelles, globular proteins 
and fibrous polysaccharides, affect a wealth of functional, 
structural and dynamical properties 0. The traditional 
approach to charged (bio)colloidal systems has been the 
mean-field Poisson-Boltzmann (PB) formalism applica- 
ble at weak surface charges, low counter-ion valency and 
high temperature [|[ . The limitations of this approach 
become practically important in highly-charged systems 
where counterion-mediated interactions between charged 
bodies start to deviate substantially from the mean-field 
accepted wisdom 0, H| . One of the fundamental recent 
advances in this field has been the systematization of 
these non-PB effects based on the notions of weak and 
strong coupling approximations. The latter approach has 
been pioneered by Rouzina and Bloomfreld Q, elabo- 
rated later by Shklovskii et al. ,7|, Levin et al. 8[, 
and brought into final form by Netz et al. [1, [H, j|. 
These two approximations allow for an explicit and ex- 
act treatment of charged systems at two disjoint limit- 
ing conditions whereas the parameter space in between 
can be analyzed only approximately [H, [l(| El, EH E3 
and is mostly accessible s olely via computer simulations 

B i M, 0, M M M E3, El • 

In the absence of a general approach that would cover 
thoroughly all the regions of the parameter space one has 
to take recourse to various partial formulations that take 
into account only this or that facet of the problem. In 
this respect the counterion-only or the one-component 
Coulomb fluid model system has proved to be of sub- 
stantial value [H . Heuristically as well as numerically. A 



proper understanding of the behavior of charged systems 
would thus start with the analysis of counter-ion distribu- 
tion around charged macromolecular surfaces, neglecting 
completely the effects of salt. 

Both the weak and the strong coupling approximations 
are based on a functional integral or field-theoretic repre- 
sentation [2(| of the grand canonical partition function of 
a system composed of fixed surface charges with interven- 
ing mobile counterions, and depend on the value of a sin- 
gle dimensionless coupling parameter 3 Q . The distance 
at which two unit charges interact with thermal energy 
U-qT is known as the Bjerrum length l^, — e\j (47T££o/cbT) 
(in water at room temperature, one has £b — 0.7nm). 
If the charge valency of the counterions is q then the 
aforementioned distance scales as o 2 £b- Similarly, the 
distance at which a counterion interacts with a macro- 
molecular surface (of surface charge density a) with an 
energy equal to k^T is called the Gouy-Chapman length, 
defined as /i = eo/(2irq£BO~). A competition between 
ion-ion and ion-surface interactions can be quantitatively 
measured with a ratio of these characteristic lengths, that 
is 3 = <7 2 £b/m = 2nq 3 £ B a/eo, which is known as the 
(Netz-Moreira) electrostatic coupling parameter [9] . The 
weak coupling (WC) regime 5 <C 1 (appropriate for low 
valency counterions and/or weakly charged surfaces), is 
characterized by the fact that the width of the counterion 
layer /i is much larger than the separation between two 
neighboring counterions in solution and thus the counte- 
rion layer behaves basically as a three-dimensional gas. 
Each counterion in this case interacts with many others 
and the collective mean-field approach of the Poisson- 
Boltzmann (PB) type is completely justified. On the 
other hand in the strong coupling (SC) regime 3 ^> 1 
(appropriate for high valency counterions and/or highly 
charged surfaces), the mean distance between counteri- 
ons, a± ~ ^/qeo/a, is much larger than the layer width 
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(i.e., a±/n ~ vS 3> 1), indicating that the counterions 
are highly localized laterally and form a strongly corre- 
lated quasi-two-dimensional layer next to a charged sur- 
face. In this case, the weak-coupling approach breaks 
down due to strong counterion-surface and counterion- 
counterion correlations. Since counterions can move al- 
most independently from the others along the direction 
perpendicular to the surface, the collective many-body 
effects that enable a mean-field description are absent, 
necessitating a complementary SC description [§[. The 
range of validity of both limiting theories at intermedi- 
ate values of the coupling parameter has been explored 
thoroughly in the literature @, H ©, [M HH, E E3 • 

Formally the weak coupling limit can be straightfor- 
wardly identified with the saddle-point approximation of 
the field theoretic representation of the grand canonical 
partition function, and is reduced to the mean-field PB 
theory in the lowest order for 5 — ► 0. The quadratic fluc- 
tuations around the mean field provide a second-order 
correction to the mean-field solution for small finite 3 < 1 
0, EJ [H, [H [3, [H, [I!]. The strong coupling approx- 
imation has no PB-like correlates Q since it is formally 
equivalent to a single particle description obtained from 
a systematic 1/2 expansion in the limit 3 — * oo, and 
corresponds to two lowest order terms in the virial ex- 
pansion of the grand canonical partition function. The 
consequences and the formalism of these two limits of the 
Coulomb fluid description have been explored widely and 
in detail (for reviews, see Refs. [!, 

Considering the inhomogeneity of charged surfaces in 
various biological contexts it has always been of interest 
to investigate not just electrostatic interactions between 
symmetrical charged surfaces, i.e. those bearing equal 
charges of the same sign, but also interactions between 
surfaces bearing unequal charges or even charges of op- 
posite sign [H, [H S H HH, HI- This problem 
has a venerable history starting from the seminal work 
of Parsegian and Gingell [13] who formulated a linearized 
PB theory of the interactions in the presence of salt. The 
linearization ansatz was later generalized in the work of 
Lau and Pincus [HI and Ben-Yaakov et al. [2!| who for- 
mulated the appropriate non-linear mean-field theory of 
non-symmetric electrostatic interactions. 

It is thus our goal in this contribution to show how 
and to what extent the asymmetry in the distribution 
of charges on two apposed planar surfaces affects the 
interactions between macromolecular surfaces carrying 
them. Below we shall present a complete analysis of 
the asymmetric case in the weak coupling limit, i.e. the 
mean-field Poisson-Boltzmann theory supplemented with 
a quadratic-fluctuations analysis, as well as in the strong 
coupling limit via the asymptotic strong-coupling the- 
ory and evaluate how these analytical results compare 
with extensive numerical simulations. We will show that 
in their respective regimes of validity (i.e. small/large 
couplings) both approximations present a very accurate 
quantitative statistical description of the system. 
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FIG. 1: (Color online) Geometry of the system comprising 
two asymmetrically charged planar surfaces located at z = ±a 
(at separation distance D = 2a) with neutralizing point-like 
counterions of valency q distributed in between [33q . 



II. GEOMETRY 

In our model system we consider uniform surface 
charge distributions on two plane-parallel surfaces (lo- 
cated at z = ±a) given by the surface charge density of 
the form 

Po(r) = o\ S(z + a) + a 2 S(z — a). (1) 
We may interchangeably use the half-separation a, or 



D = 2a 



(2) 



to identify the surface-surface distance. 

We assume furthermore that the charge of both bound- 
ing surfaces is compensated by mobile counterions of 
charge valency q immersed in an aqueous medium of di- 
electric constant e and distributed in between the two 
surfaces (see Fig. [1]). We thus neglect all coions. This ap- 
proximation is relevant for low salt concentrations where 
the Debye screening length is much larger than the scales 
of interest [34[ . We consider the surfaces as impenetra- 
ble to counterions and neglect the dielectric discontinuity 
across the bounding surfaces which was addressed at var- 
ious levels of approximation in [H, [H, H3, HI] . 

Without loss of generality we can assume here that 
q > and 

°i + °2 < 0, and (72 > o"i, so that <j\ < 0. (3) 

It will be helpful for our later developments to introduce 
an asymmetry parameter £ that will allow us to quantify 
the dissimilarity between the two bounding surfaces as 



-1. 



(4) 



Furthermore, by suitably normalizing the results one can 
concentrate exclusively on the interval — 1 < £ < 1. All 
other cases can be mapped onto this interval with appro- 
priate rescaling of the parameters. The values £ = 1 and 
£ = — 1 represent exceptional points in the parameter 
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space: £ = 1 is the standard symmetric case (<ti = 02) 
already amply treated in the literature, and £ = — 1 rep- 
resents the antisymmetric case (<ji = — a 2) with no coun- 
terions between surfaces that reduces to the trivial case of 
a planar capacitor. These two well-understood limiting 
cases will be thus omitted from our discussion. 

Counterions between the surfaces satisfy electroneu- 
trality condition that can be written in the form 

Ne q+(a 1 +a 2 )S = 0, (5) 

where S denotes the (infinite) area of each surface. 

III. DIMENSIONLESS REPRESENTATION 

Because of the asymmetry present in the system, we 
have two length scales describing the interaction of the 
counterions with each of the bounding surfaces. These 
two length scales are given by the corresponding Gouy- 
Chapman lengths associated with the two surfaces as 



/'1 



M2 



2n£ B q\cT2\ C 



(6) 



For the same reason we can thus define two different cou- 
pling parameters 



g 2 e B _ p 

Mi 



^2 



(7) 



each one being defined by the ratio between the Bjerrum 
length and the corresponding Gouy-Chapman length. In 
what follows, we rescale the surface separation as 



D = D/n 



(8) 



(or the rescaled half-distance as a = a/fi) with respect 
to plate 1. With an appropriate rescaling one could also 
equivalently define all dimensionless lengths with respect 
to plate 2. Thus the minimal set of dimensionless param- 
eters that fully characterize the system in the thermody- 
namic limit is given by {S, D}. 

Other physical quantities such as the mean electro- 
static potential ^p(z), number density of counterions, 
n(z), and the pressure, p, acting on each surface can be 
rescaled as well. We shall use the standard rescaled elec- 
trostatic potential 

as well as the rescaled density and pressure 
n(z) j (3p 



n(z) 



27r€ B (CTi/e ) 2 



and p 



(9) 



where f3 = 1/ k B T, and all the other quantities have been 
defined above. 



IV. MEAN-FIELD POISSON-BOLTZMANN (PB) 
APPROXIMATION 



In the weak coupling regime, the leading contribution 
to the partition function comes from the saddle-point 
configuration of the local fluctuating electrostatic poten- 
tial, tpo(z) 20}. The saddle-point configuration can be 
straightforwardly translated into a solution of the PB 
equation and corresponds to an exact asymptotic result 
in the limit S — > [a, |23[ . For a system containing only 
counterions, the PB equation for the dimensionlesspo- 
tential, ^>q(z), can be written in the standard form 



d 2 Vo(z) 



dz 2 

with boundary conditions 

dtpp 
dz 

d^o 
dz 



(4ir£ B q 2 )\ e-^ z \ 



(10) 



2 



A 4 ' 



(11) 



Integration of the PB equation gives rise to the first in- 
tegral of the system of the form 



(3po 



1 (diPo 
'8n£ B q 2 \ dz 



n Q (z), 



(12) 



where the constant po is nothing but the mean-field PB 
pressure acting between the bounding surfaces and 



n (z) = \ e~^ z \ 



(13) 



is the PB number density profile of counterions between 
the surfaces. The normalization factor Ao follows from 
the electroneutrality condition ([5]) as 



A n = - 



0\ + (T2 



qe J° a dze-+»( z ) 



(14) 



The nature of the solution ipo(z) obviously crucially 
depends on the sign of the pressure po [H, [2i| . Different 
forms are obtained for positive and negative pressures, 
corresponding to repulsion and attraction between the 
bounding surfaces respectively. We review these different 
cases separately. 



A. Repulsion regime po > 

In the case of repulsive pressure the appropriate solu- 
tion of Eq. pop can be written as 



tpo = In I cos 2 a{z- z )\. 



(15) 
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where the constants zq and a 2 = 27r^B<3 ,2 (/3po) are ob- 
tained from the boundary conditions and satisfy the 
set of two equations 



atano(a + Zq) = 



1 



atana(a - Zq) = — 



(16) 
(17) 



Eliminating z we obtain an equation for a of the form 



tan (2aa) 



ory? — C 



(18) 



The solution of this equation provides the final result for 
the repulsive PB pressure pq. Note in particular that in 
rescaled units and by definition one has 



Pa = a 2 , 



(19) 



where a — afi. Once a is known, the parameter zq can be 
simply obtained from Eqs. (fTrj|) or (JT7J) and thus the po- 
tential ipo(z), Eq. (fTS")) , is fully determined. The density 
profile of counterions then follows from 



n (z) 



0Po 



cos 2 a(z — Zq) 



(20) 



A positive (repulsive) solution for the pressure as con- 
sidered in this section is always possible for any given 
asymmetry parameter £ (excluding the trivial case of 
£ = —1). In particular, it easily follows that within the 
mean-field theory two surfaces of equal sign (£ > 0) al- 
ways repel, that is at all separation distances D. When 
the surfaces bear charges of opposite sign £ < 0, they 
attract at large separations (see below) and a repulsion 
emerges only at sufficiently small separations. 

At small separations D -C 1, we can obtain the lim- 
iting solution for a and thus the limiting small-distance 
pressure po as 



Po(D) 



i + C 
D 



(21) 



for arbitrary |£| < 1 as noted in Scetion [TTJ This is 
of course nothing but the ideal-gas osmotic pressure of 
counterion confined between the two plates (i.e., po = 
Nk-B,T/(SD) in actual units), which dominates over the 
energetic contributions at small separations. At large 
separations D » 1 and for £ > 0, we obtain the asymp- 
totic expansion 



Po(D) = w 



(D 



\D 2 



(22) 



which is valid for D > 2(1 + l/£). Thus in the limit 
D — ► oo, the pressure behaves as 



Po(D) c w 



(23) 



which agrees with the asymptotic pressure between two 
equally charged surfaces (£ = 1). For smaller £ than 
determined above, i.e. for 3> D » 1, one needs to 
invoke a different asymptotic expansion and specifically 
for £ ~ (one surface being neutral), one obtains 



Po(D) 



AD 2 



(24) 



This latter asymptotic result may be obtained from the 
one in Eq. ([23| by redefining D — > 2D. This may be 
understood simply by noting that because of symmetry 
a system with £ > may be decomposed into two halves 
each with an effective asymmetry parameter £ = 0. 



B. Attraction regime po < 

An attractive pressure on the mean-field level is pos- 
sible only if the surfaces are oppositely charged £ < 0. 
The appropriate solution in this case is given by 



'•'(! = In<j -^-j- sinh 2 a(x - ;,,)}>. 



(25) 



where the constants zq and a 2 = 27t£b9 2 (/3|po|) can again 
be obtained from boundary conditions, this time in the 
form 



acotha(a + z ) 
a coth a(a — zq) 



1 

H 

C 



Eliminating zq we obtain an equation for a as 

coth (2aa = -±-JL—. 

txa(l + C) 

In this case we have in rescaled units 



po = -a 2 , 



and for the density profile of counterions 

0\po\ 



n (z) 



sinh a(z — Zq) 



(26) 
(27) 

(28) 
(29) 
(30) 



The asymptotic form of the attractive pressure at large 
separations D ^> 1 can be derived as 



p Q (D)^-C U-4 



1 + C 

i-C 



,2(D 



(31) 



where ( < as noted above. For infinite separations, this 
pressure does not vanish and exponentially approaches 
— C 2 since for £ < the system behaves partially as a 
simple capacitor. 
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C. Zero pressure po = 

In the case of charged surfaces with opposite sign 
(C < 0), the large-distance attraction regime and the 
short-distance repulsion regime merge at the point of zero 
pressure, D — D* , where the surfaces are at equilibrium. 
In this case, the PB solution for the potential reads 



■ipo = ln|27r^ B g 2 A (z 



zo) 



and the density profile of counterions is given by 

1 



n (z) 



2n£ B q 2 {z - z ) 2 ' 



(32) 



(33) 



where zo is found from the boundary conditions as zo = 
—fi — D* /2, and the bound-state separation, D* , follows 
in rescaled units as 



D* = - 



i + C 

c ■ 



(34) 



The surfaces attract for D > D* and repel for D < D* . 
In the vicinity of D*, that is for \D* — D\ < 1, the 
pressure behaves as 



Po 



3C 3 



i + C 3 



\D* - D\ 



(35) 



This concludes the calculation of the inter-surface pres- 
sure on the mean- field PB level strictly valid for S — > 0. 

The preceding results may be summarized in a phase 
diagram shown in Fig. 3] in terms of D* and the asymme- 
try parameter £ displaying the mean-field attraction and 
repulsion regimes separated by the boundary line (|34[) . 
The forms of the pressure here are completely consistent 
with those derived by Lau and Pincus [28j via a different 
route. 



V. WEAK-COUPLING (WC) ANALYSIS: 
QUADRATIC FLUCTUATIONS AROUND MEAN 
FIELD 

The first non-zero correction to the saddle point is sec- 
ond order in the fluctuations of the local electrostatic po- 
tential around the mean-field PB solution, i/V Our goal 
here is to calculate the corrections in pressure, p2(D), 
stemming from these quadratic fluctuations, which leads 
then to the total WC pressure 



p(D) = p (D)+p 2 (D). 



(36) 



This approach has correlates in many diverse areas of 
physics where fluctuations around a mean-field solu- 
tion are important [39| and goes under different names, 
though the physics is always the same. We may conven- 
tionally refer to the mean-field PB term, p$(D), and the 
fluctuations contribution, p2(D), as the zeroth- order and 



the second-order correction terms on the WC level, re- 
spectively. This procedure formally also corresponds to 
a series expansion in powers of S (loop expansion) around 
the asymptotic mean- field solution (S -> 0) [1 So, [II, [1| 
and is thus expected to be valid for sufficiently small cou- 
pling parameters as will be determined later. Note also 
that in this latter sense the second-order pressure turns 
out to be proportional to S, that is m ~ S, and thus 
corresponds to a first-loop correction [i| [23| . 

In order to proceed, one needs to evaluate the appro- 
priate Hessian of the field action in the partition function 
and study its fluctuation spectrum (see Refs. [2(| H(| for 
more details). The Hessian of the field action can be 
derived in the form 

H(r, r') = u-^r, r') + /3(e 0< z) 2 n (z) S 3 (r - r'), (37) 

where w _1 (r,r') = — ££oV 2 (5 3 (r — r') is the inverse 
Coulomb operator and Hq (z) is the zeroth-order PB den- 
sity as derived in the previous section. Hence, 



(4ttW)M 2 ) 



2a 2 



cos 2 a(z — Zo) 



(z - z ) 2 
2a 2 



sinh a(z — zo) 



Po > 0, 
Po = 0, 
Po < 0. 



(38) 



The corresponding correction, T2, to the free energy of 
the system is then given by the trace-log of the Hessian. 
It can be written equivalently in the following form (20l . 

/W-ilk h *(,.0-!jf Ob^dO. (38) 

This form can be derived rather straightforwardly by us- 
ing the argument principle (2~0 | and converting the dis- 
crete sum of eigenvalues of the Hessian operator into an 
integral over the transverse wave-vector Q = {Q x ,Qy), 
with density of modes S/ (2-7r) 2 of the logarithm of the 
secular determinant T>\ of the same operator. The in- 
dex A in the secular determinant refers to the eigenvalue 
equation that can be derived in the form 

- Q 2 - A(4^b<Z 2 HW)/a(Q, *) = 0. (40) 

By simply writing f(a, Q) for the quotient 
T>i(Q)/T>o(Q), and noting that the secular deter- 
minant depends explicitly also on the value of the 
inter-surface spacing, a, the free energy contribution 
from the quadratic fluctuations can be equivalently 
expressed exactly in a dimensionless form as 



T2 
S 



QlnV(a,Q)dQ. 



(41) 
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where T% — 2eeoJ r 2/crffJ, 3 , Q = [iQ and the rescaled area 
S = S/n 2 . Here the secular determinant of the Hessian 
for homogeneous transverse modes has been written as a 
function of dimensionless quantities a — a/fi, Q = [iQ. 
This determinant has to be standardly regularized so that 
all irrelevant constants, i.e. all the terms not depend- 
ing on the separation between the bounding surfaces, are 
dropped, amounting to a rescaling 



V(a,Q) 



T>(a — > oo, Q) 



(42) 



This corresponds to a subtraction of the part of the free 
energy for two separate interfaces at infinite separation 
from the total free energy. 

In the next step one has to calculate the secular de- 
terminant V{a,Q) for each of the pressure regimes sep- 
arately, since the appropriate cigenfunctions of the Hes- 
sian depend on the mean-field solution that in its turn 
depends on the sign of the interaction pressure, see Eq. 
([38| . In what follows we shall follow closely the deriva- 
tions in Refs. HE El Si]- 

The total pressure in the weak-coupling limit is thus 
the sum of the PB pressure and the quadratic fluctuations 
correction and can be written as 

p(D) = p (D) +P2(D) = Po(D) - i fe J • (43) 



A. Repulsion regime po > 

In this regime the secular determinant of the Hessian, 
Eq. (|4H|) . can be obtained by solving 



dz 2 



2a 2 



cos 2 a(z — zq) 



y(Q,z) = 0, (44) 



with appropriate boundary conditions implying continu- 
ity of the solution and its derivative across the bounding 
surfaces at z = ±o. 

The general solution of Eq. |44|) for various regions in 
the perpendicular direction can be written in the form 



y{Q,z) 




(45) 



where 



yi 

2/2 



i X 

1 + — tana(z — Zq) 



-Qz 



Q 2 



l 



Q 



tana(z — z ) 



(46) 



(47) 



Taking into account the continuity of the solution and its 
derivatives we get a set of four homogeneous equations for 



the coefficients A, B, C and D. The solution exists only 
if the (secular) determinant of this system equals zero. 
Thus we derive the secular determinant of the Hessian 
operator in this case in the rescaled form 



V(a,Q) 



r+(fi,Q)-r+(a,0)e- 



where 



T + (a, Q) = {l+& 2 + 2Q + 2Q 2 )(( 2 



a 2 



-4Q5 



2(Q 



(48) 



2Q 2 ). 

The regularized form of the secular determinant is ob- 
tained by taking the quotient as indicated in Eq. (|42|). 
While doing this, it is important to realize that a also 
depends on the inter-surface distance. In fact from Eq. 
(f2"3")l it follows that the appropriate limit of a is 



Jim a(a) = 0. 



(50) 



In the regularization of the secular determinant this lim- 
iting behavior should be consistently taken into account. 

Finally, the dimensionless quadratic fluctuations free 
energy T% can be calculated numerically via Eqs. (|41[) 
and (|48|) . The fluctuations contribution to the pressure, 
p~2, then follows from Eq. (|4"3")l . 

The asymptotic form of the second-order dimensionless 
pressure can be obtained analytically. Note that at large 
separations D — > oo, a repulsive mean-field pressure po > 
0, as considered in this section, is possible only for non- 
negative £• For not too small £ > 0, i.e. when D ^S> 
2(1 + 1/C), we find 



p 2 (D) 



-Zir 2 



hxD 



while for £ ~ 0, i.e. when ( 1 3> D ^> 1, we get 



P2(D) 



-Zir 2 



InZ) 



(51) 



(52) 



Again the difference in the two cases above is due to the 
symmetry of the problem in the latter case, that can be 
described by redefining D — > 2D and discarding the sub- 
dominant terms. 

The second-order pressure is obviously attractive and 
in this regime leads to a reduction of the total pressure 
from the mean-field value po- This clearly shows that 
electrostatic correlations favor attraction between two re- 
pelling asymmetrically charged plates. However, the to- 
tal pressure never becomes negative as the fluctuations 
are assumed to be small within the second-order weak- 
coupling analysis. 



B. Attraction regime po < 

In this case the secular determinant of the Hessian, Eq. 
P0|) . is obtained by solving 



2a 2 



dz 



sinh ol(z — zq) 



y(Q,z) = Q. (53) 
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The general solutions for particular regions in the z di- 
rection are 



y(Q,z) 




(54) 



where 



yi 



1)2 



1 — — cotha(z — Zq) 



,-Qz 



1 + — cotha(z — zq) 



(55) 
(56) 



Again the solution exists only if the determinant of 
the system of equations, which connect coefficients 
E, F, G, H and stems from the application of the bound- 
ary conditions at z = ±<z, is identically zero. This again 
defines the secular determinant T>(a, Q) appropriate for 
this case. It is easy to show that the secular determinant 
T>(a, Q) can be obtained from the po > result, Eq. (|4"8")l . 
simply by substituting a 2 — > —a 2 and so we can write 



r_(<s,Q) -r_(a,o)e- 4<5a 



Q 2 -& 



■ \ — 



(57) 



where 



r_(a, Q) = (1 - a 2 + 2Q + 2Q 2 )(( 2 - a 2 + 2(Q + 2Q 2 ). 

(58) 

Here we can again regularize the secular determinant to 
discard divergences, Eq. (14"2")) . Again one has to be care- 
ful by taking the correct limit for a in the above regular- 
ization scheme. In this regime, the appropriate limit is 
given by 



Jim <5(a) = — £, 



(59) 



as follows straightforwardly from Eq. (|31j) . The fluc- 
tuations contribution to the pressure, p2, can then be 
evaluated numerically from Eq. (|43[) . 

The asymptotic form of f>2 for D ^S> 1 can be derived 
analytically as 



p 2 (D)~Zf(Qe 



2C.D 



(60) 



which is applicable only for charged surfaces of opposite 
sign, £ < 0, which can attract (po < 0) at large separa- 
tions. The function /(£) is defined as 



/(C) = c 



3 1 + C / 2 arctan y/l - 2( 2 




In 



>3 1 + C / 2 tanh -1 v/2( 2 - 1 



i-C 



111 



i-C 2 

2C 2 



i-C 2 

2C 2 



(61) 



(62) 



for — 1 < ( < — V2/2. The second-order pressure thus 
asymptotically decays exponentially and can be only at- 
tractive. For not too small values of C, it is thus qualita- 
tively very different from the case £ > 0. 

The total weak-coupling pressure, p = P0+P2, is shown 
in Fig. [5] for a few different asymmetry parameters and 
a relatively small value of the coupling parameter S. 




FIG. 2: (Color online) Rescaled weak-coupling inter-surface 
pressure, Eqs. (|19|l and (|29[) , as a function of the rescaled 
distance, D, between two charged plates for three different 
asymmetry parameters £ = —0.5, and 0.5 as shown on the 
graph. The total pressure, p = po + f>2 (dashed lines, plotted 
here for S = 1 using Eq. (J43J ) , is always lowered from its 
mean-field PB value, po (solid lines, obtained for H — > 0), since 
quadratic fluctuations around mean field favor attraction. 



C. Regime of validity of the weak-coupling theory 

As noted above the foregoing weak-coupling analysis is 
valid as long as the quadratic corrections are sufficiently 
small so that the series expansion around the mean-field 
solution does not diverge [§, [23|. As an approximate 
measure for the validity regime of this scheme, one can 
require that the second-order correction term is smaller 
than the leading order term, i.e. 



(63) 



This leads to a useful criterion identifying the regime of 
coupling parameters and distances in which the weak- 
coupling theory is applicable. For po > and by employ- 
ing the closed-form expressions obtained for large sepa- 
rations D 3> 1, we find the validity condition 



a < 



D 
ItZd' 



(64) 



This indicates that at a given non- vanishing 5, the weak- 
coupling scheme becomes increasingly more accurate at 
larger separations, while as the surfaces get closer a 
smaller coupling parameter needs to be chosen. 
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On the other hand, for po < (which occurs for £ < 0) 
and at large separations D 1, we obtain 



a < 







1/(01 



-2C-D 



(65) 



The right hand side here is exponentially large mean- 
ing that for charged surfaces of opposite sign, the weak- 
coupling analysis performs far better at finite coupling 
parameters and smaller inter-surface separations than for 
the surfaces of equal sign (£ > 0). 

Finally, note that for po = that corresponds to the 
equilibrium phase boundary line in Fig. |4] for ( < 0, 
we deal with a situation where the leading order term 
is zero and the fluctuations are dominant at any finite 
value of 2. The convergence of the loop expansion has to 
be determined in this case by evaluating the higher order 
terms which we shall not consider in this paper. 



Differentiating the free energy with respect to the 
surface-surface distance D we get the corresponding pres- 
sure acting between the bounding surfaces 



p(£) = -I(l + C 2 ) + i(l-C a )coth 



(W) : 



D 



(71) 



The dependence of this dimcnsionlcss pressure on the 
separation for different values of ( is presented in Fig. 
[31 Note that the SC pressure can become attractive for 
both like-charged and oppositely charged surfaces which 
contrasts with the mean-field theory that does not allow 
attraction between like-charged surfaces. This is because 
of the strong electrostatic correlations mediated by coun- 
terions between the charged surfaces for 3 ^ 1 and has 
been investigated throughly before for equally charged 
surfaces @. Our results show that a similar attraction 
mechanism holds for asymmetrically charged surfaces in 
the SC limit. 



VI. STRONG-COUPLING (SC) THEORY 

The strong-coupling approximation coincides with the 
lowest order non-trivial expansion of the partition func- 
tion in terms of the fugacities of the counterions. This 
expansion may be expressed as a 1/2 series expansion Q, 
whose leading order term (2 — * oo) corresponds to the 
so-called SC theory. We will not delve into the strong- 
coupling expansion in more detail since it has been ex- 
haustively reviewed in the literature 0, IE 0- On the 
leading order, the free energy is obtained as 

T = W - Nk B T\n [ e -^ Wl+w ^dV, (66) 




where Wo is electrostatic interaction energy of charged 
surfaces 



Wo 



2ee 



SD, 



(67) 



with S representing surface area of each plate, and W\ 
and W2 are electrostatic interaction energies between a 
single counterion and individual charged surfaces, i.e. 



*~^<f + .). 



w 2 



qe Q a 2 ,D , 



(68) 



Since in the strong-coupling regime the free energy is 
given via simple quadratures, it is much simpler to eval- 
uate it than on the weak-coupling level. Defining the 
rescaled free energy 



- 2eeo ^ 
T = -T 



we obtain 



J = (l + C 2 )f -(1 + In sinh 



u-o? 



(69) 



(70) 



FIG. 3: (Color online) Rescaled strong-coupling inter-surface 
pressure, Eq. (|71[) , as a function of the rescaled distance, D, 
between two charged plates for different asymmetry param- 
eters £ = —0.5, 0, 0.5 and 1 as shown on the graph. Long- 
distance attractive pressure is suppressed for the £ = case. 
For symmetrically charged surface (£ — 1), we recover the 
standard SC result p(D) = -1 + 2/.D 0. 

The pressure exhibits two well-defined limiting laws 
obtainable for small and large inter-surface separations. 
For small separations D <C 1, we have 



p(D) 



1 + C 
D ' 



(72) 



for arbitrary |£| < 1 as noted in Section [TTJ For large 
separations D 3> 1, we obtain 



p(D) ~ -C 2 



(73) 



Note that in the limit D — > the SC pressure coincides 
with the PB result on the leading order (compare Eqs. 
(|2~Tj) and l[72p) and represents the ideal-gas osmotic pres- 
sure of counterions which dominates over the electrostatic 
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contributions. The PB and SC forms for £ < coincide 
also in the limit of D ^> 1 (compare Eqs. (f5Tj) and ([75]) ) 
and reduce to the pressure in the capacitor. 

The dependence of the pressure on the inter-surface 
separation points to the existence of a bound state de- 
fined via p(D*) = 0. The SC bound-state separation D* 
can be expressed analytically as 

21n|C| 



D* 



i-C 



(74) 



and is presented in Fig. U as a function of £. Obvi- 
ously both like-charged and oppositely charged surfaces 
can form bound states at small surface-surface separa- 
tions. The bound-state separation approaches infinity 
and the surfaces unbind asymptotically as £ — > 0, that is 
when one plate becomes electroneutral. 



PB 

WC (2nd, S = 1) 
SC 



' \ SC attraction 




v WC repulsion 



FIG. 4: (Color online) Rescaled bound-state surface- surface 
separation D* = 2a* as a function of the asymmetry pa- 
rameter £ as predicted by the zeroth-order PB theory (solid 
line), Eq. (|34jl . the second-order WC theory at 3 = 1 
(dotted line), and the strong-coupling theory (dashed line), 
Eq. (|74[) . Surfaces attract for D > D* and repel other- 
wise. Only in the SC limit can two surfaces of equal sign 
(( > 0) attract. The second-order WC result is obtained 
by finding the zero-pressure point of the total WC pressure 
p(D) = po(D) +p 2 {D), Eq. |g3), for given 3 and C- 



as a function of z and £. 



A. Regime of validity of the strong-coupling theory 

The regime of applicability of the leading order SC 
theory follows from a simple criterion that has been dis- 
cussed and confirmed previously in the case of equally 
charged surfaces by both MC simulations and higher- 
order calculations [3, [t| . The generalization to asym- 
metrically charged surfaces is straightforward. 

For large couplings, counterions are strongly attracted 
to an oppositely charged surface as the counterion-surface 
interaction becomes large and equivalently, the Gouy- 
Chapman length, fi, and the thickness of the counteri- 
onic layer at the surface become small. The layer thick- 
ness has to be compared with the typical lateral spacing 
between counterions aj_ ■ For counterions sandwiched be- 
tween two asymmetrically charged surfaces, this latter 
quantity follows from the local electroneutrality condi- 
tion as 



<?e 



(77) 



up to a factor of the order unity and assuming that the 
surfaces are sufficiently close so that they may be strongly 
coupled via the counterions as will be determined consis- 
tently here. In rescaled units, one gets 



a 2 , 



i + C 



(78) 



where a± = a±/fi. Obviously, a± becomes large relative 
to the layer thickness as S grows. Note that on the other 
hand the lateral Coulomb repulsion between counterions 
in this quasi-two-dimensional layer becomes much larger 
than the thermal energy, i.e. q 2 iB/a± ~ \/S ^> 1, indi- 
cating that counterions form a strongly correlated liquid 
in which they are highly localized within correlation holes 
of lateral size aj_ ~ VE ^> 1 0, la, Q . Thus for surface- 
surface separations, D, smaller than the correlation hole 
size, i.e. 



Finally one can also derive the explicit form of the 
counterion density as a function of the normal coordi- 
nate z. This can be read off simply from the integrand in 
Eq. ([66]), that is n(z) = Cexn(-0(Wi + W 2 )), where 
C is a normalization factor According to the 

electroneutrality condition, we normalize the density as 
J" n{z) dz = — (<7i + o^Vqeo, or in rescaled units 



h{z) dz = 1 + £, 



(75) 



where we have defined z = zj /i with fi — ~e$/ (2tt£bQO'i) 
being the Gouy-Chapman length with respect to plate 1. 
From here, the density profile is obtained as 

i - C 2 e -( 1 -o^ 

= — sinh [(1 - C)D/2] ' (?6) 



i + C 



(79) 



counterions can move almost independently from each 
other in the direction normal to the surface and one can 
safely assume that the effective surface-surface interac- 
tion as well as the counterionic density profile follow only 
from the interactions of individual counterions with the 
bounding charged surfaces. The counterion-counterion 
interactions contribute on the sub-leading order and mat- 
ter at larger separations. This picture is of course con- 
firmed on a systematic level by the SC expansion anal- 
ysis d, H, Q , and the above equation sets a criterion for 
the validity regime of the single-particle leading order SC 
theory (3 — > oo), when applied to finite coupling param- 
eters. 
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VII. SIMULATIONS 

We performed Monte-Carlo (MC) simulations in order 
to study the system of two asymmetrically charged sur- 
faces beyond the analytical limits of weak and strong cou- 
pling discussed above. All simulations were performed 
in the Canonical ensemble (NVT) using the standard 
Metropolis algorithm [Z(|. The mobile counterions were 
modeled as point charges [33| enclosed in a simulation 
box bounded in z direction by two charged surfaces of dis- 
tance D and surface charge densities a\ and 01 (compare 
Fig. QJ. Periodic boundary conditions were applied in the 
lateral directions parallel to the bounding surfaces. The 
lateral size of the charged surfaces, L, which is set equal 
to the lateral size of the simulation box, was held fixed 
throughout all the simulations. The number of counteri- 
ons, N, was varied between 60 and 1800, depending on 
the system parameters in order to fulfill the electroneu- 
trality condition. The counterions interact through the 
Coulombic potential 



U{n 3 ) 



(80) 



with qi = e^q being the charge of the i-th counterion 
and Tij the separation distance between the i-th and the 
j-th counterions (with i, j = 1, . . . , N). The interaction 
energy of the i-th counterion with the charged surfaces 
in the simulation box are given by 



U k {zik 



AL In 



arcsin 



VP? 
(^72)4 



[(W+41 2 



7T 

2 



where — o\ and oi is surface charge density of the k- 
th surface (with k = 1, 2) and Zik is the normal distance 
between the z-th counterion and the A:-th surface. The 
long-ranged Coulomb interactions in this system were ac- 
counted for via a charged sheet scheme similar to that 
proposed by Torrie and Valleau [4lT | . This scheme makes 
use of the counterion profile in the simulation box in order 
to calculate an external field, stemming from the long- 
ranged interactions. This external field is iteratively up- 
dated and self-consistency is achieved normally in a few 
iterations. 

In the course of simulations, new configurations were 
created by trial displacements of the counterions and 
equilibration was accomplished by running through 10 6 
configurations. The 10 7 following configurations were 
then used for the production runs. 

The pressure p was calculated in the production runs 
according to the contact-value theorem as 



Pk = k B Tn c r tact 



2ee ' 



(82) 



where n™ ntact is the density of counterions at contact 
with the fc-th surface. In thermodynamic equilibrium, 



the pressure does not depend on which surface [k = 1 or 
k = 2) is chosen in order to calculate the pressure from 
the above equation, and the contact condition at both 
surfaces leads to precisely the same value for the pres- 
sure. All simulations were conducted at fixed tempera- 
ture T — 298 K, lateral simulation box size L = 245 A 
and dielectric constant e — 78.7, which is assumed to be 
the same throughout the system. 

The simulations were performed at different values of 
the coupling parameter 5 and the asymmetry parameter 
£. We explored the 5 parameter space extensively by us- 
ing S = 0.32, 0.64, 2.5, 3.2, 5.1, 6.4, 8.6, 17, 25, 51, 86 and 
172 in order to cover exhaustively both the weak cou- 
pling and the strong coupling regimes. The concurrent 
values of the asymmetry parameter were always taken 
as C = —0.5, 0, +0.5 at each value of the coupling pa- 
rameter. The results are plotted in the form of rescaled 
density and pressure as previously defined in this paper 
(Section [m|. 



VIII. DISCUSSION 

In order to asses the validity of the weak and strong 
coupling results presented above for asymmetrically 
charged surfaces, we performed extensive MC simula- 
tions and compared them to analytical results in both 
limits. It transpires from this comparison that the sim- 
ulation results corresponding to an exact evaluation of 
the partition function are always bracketed by the WC 
and the SC limiting forms, smoothly approaching them 
in the appropriate limits of the coupling parameter S. 

First we compare the density profiles of simulations 
with theoretical results given by Eqs. (|20p and (|30| for 
the PB limit (solid lines) and Eq. ([75]) for the SC limit 
(dashed lines). As seen from Fig. [5] the theoretical PB 
and SC rescaled density profiles represent two extremal 
cases and all MC simulations results with finite values of 
5 are located consistently between these two limits. The 
MC results for small S are almost exactly spot on the PB 
prediction, while larger discrepancies are observed as S 
grows. For large enough 5 > 10, the MC results slowly 
converge to the SC result. This is especially clear for 
surfaces with charges of equal sign, £ = 0.5, whereas for 
surfaces with opposite sign, £ = —0.5, there is no big 
difference between PB and SC profiles. 

Next we consider the inter-surface pressure as obtained 
from the simulations (symbols in Fig. [BJ as well as the 
PB theory, Eqs. (JTSJ and [[25]). and the SC theory, Eq. 
(jTTj) (solid and dashed lines, respectively). The PB re- 
sult is expected to be valid for separations i) > H || (see 
also Section ESI)- Therefore, for 2 = 0.32 the PB line 
expectedly agrees nicely with the simulation data (open 
squares) in the whole range of separations shown in the 
figure. Upon closer inspection, however, we find small 
deviations from the PB result as shown in the insets in 
Fig. [5]for all three values of the asymmetry parameter £. 
In this case, the fluctuation correction to the mean-field 
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pressure accurately compensates for these deviations and 
adding the second-order correction to the PB pressure 
leads to a total WC pressure, Eq. (|43|) (shown as a dot- 
ted line in the inset) that matches the simulation data 
perfectly. Note that the pressure changes can be dras- 
tic even on the WC level as C assumes different values. 
For £ = 0.5 and the WC pressure is strictly repulsive, 
while it turns attractive and leads to a bound state (zero 
pressure point) for £ = —0.5. 

In the intermediate regime of coupling parameters, the 
simulation results for the pressure are clearly bracketed 
by the two limiting analytical forms, given by the PB plus 
the second order correction and the SC expressions of the 
interaction pressure. The SC prediction is expected to 
be valid for separations D <C vH as discussed in Sec- 
tion IVI Al Consistently, the interaction pressure starts 
off close to the strong coupling limit at small separations 
and then smoothly converges to the weak coupling limit 
for larger separations. This is strictly true for Q = 0.5 
and 0. In the case of £ = —0.5 the difference between 
the strong and weak coupling results for the rescaled in- 
teraction pressure is marginal and the simulation data 
and the analytical results nearly coincide for all rescaled 
separations D = D//J,. We emphasize that the pressures 
and the density profiles are plotted here in rescaled repre- 
sentation; in actual units, Fig. [^corresponds to different 
ranges of separation, D, for the WC and SC regimes as 
the Gouy-Chapman length, /i, is typically very different 
between the two limits (small at high couplings and large 
at small couplings as may be realized, e.g., by changing 
the counterion valency at fixed surface charge densities 
and Bjerrum length). 

For large values of the coupling parameter the sim- 
ulation results for the interaction pressure expectedly 
follow very closely the strong coupling prediction for a 
wider range of inter-surface separations. The correspon- 
dence between the SC theory and simulations is better 
for £ = —0.5 than for ( = 0.5, which can be again traced 
back to the fact that the strong and the weak coupling re- 
sults are very close to one another for charged surfaces of 
opposite sign in the whole range of rescaled separations, 
whereas they differ significantly in the case of surfaces of 
equal sign. For intermediate and large couplings, we have 
not attempted to compare our data with the second-order 
WC approximation as this approximation breaks down at 
the range of distances shown in the figures (Section lV C| . 

Note also that in all cases considered here the theoret- 
ical and the simulated values of the interaction pressure 
converge for very small surface-surface separations. In 
fact, the SC and PB results coincide in the leading or- 
der as the rescaled distance, D, tends to zero, Eqs. (j2"Tj) 
and |72|) . since both are dominated by the osmotic pres- 
sure of counterions. The sub-leading corrections for very 
small D are different in the PB and SC limits and on 
this level the simulation data with finite 3 are generally 
expected to agree better with the SC prediction at small 
separations @. 

We now consider the simulated bound-state separa- 
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FIG. 5: (Color online) Rescaled counterion density profile 
h(z) between two asymmetrically charged surfaces at half- 
separation 5 = D/2 = 1.34 for three different asymmetry 
parameters C, = 0.5,0 and —0.5 (top to bottom). Solid lines 
represent the PB prediction, Eqs. (|20|l and (|30|l . and the 
dashed lines show the SC prediction, Eq. (|76|l . Symbols cor- 
respond to MC simulations data at three different coupling 
parameters H = 0.32 (open squares), 8.6 (filled squares) and 
86 (open circles). 



tion, D* , as a function of the coupling parameter in Fig. 
[7] (symbols). For £ — 0.5, the PB theory (S — + 0) gives 
only repulsion and predicts no bound state. While the 
SC theory predicts a closely packed bound state with D* 
given by Eq. (fTI)) . explicitly here D* ~ 2.77. As seen, 
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FIG. 6: (Color online) Rescaled interaction pressure, p, as a 
function of the rescaled inter-surface distance, D, for three 
different values of the asymmetry parameter £ = 0.5, and 
—0.5 (top to bottom). Solid lines represent the PB predic- 
tion, Eqs. (|19[) an d l|29[) , dotted lines show the second-order 
weak-coupling pressure (PB plus second-order corrections), 
Eq. (|43[) . an< A t ne dashed lines are the SC prediction, Eq. 
(|71|) . Symbols correspond to MC simulations data at three 
different coupling parameters S = 0.32 (open squares), 8.6 
(filled squares) and 86 (open circles). Insets show details at 
small pressures along with the PB theory result as well as the 
second-order WC result for H = 0.32. 



by increasing the coupling parameter the simulation data 
for D* (filled squares) decrease monotonically and rather 
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FIG. 7: (Color online) Rescaled bound-state separation D* 
as a function of the coupling parameter S for different values 
of f. Main set: symbols are simulation data for £ = 0.5 (filled 
squares) and —0.5 (filled circles). The SC results, Eq. (|74)| . 
for £ = —0.5 and 0.5 are represented by dashed and dotted 
lines, respectively. The inset represents the detailed view for 
£ = —0.5, where we also show the zeroth-order PB result, Eq. 
(|34[) (solid line), and the result from the second-order WC 
approximation (dot-dashed line) obtained numerically from 
Eq. J43). 
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FIG. 8: (Color online) Rescaled bound-state separation D* 
as a function of £. Symbols are simulation data for different 
values of the coupling parameter at £ = 0.5 and —0.5. The 
region around £ = —0.5 is expanded in the inset showing the 
crossover from the PB prediction (solid line) to the SC pre- 
diction (dashed line) upon increasing the coupling parameter. 



slowly converge to the SC prediction (dotted line). The 
comparison is again worse for the £ = 0.5 than for the 
C = —0.5 case (dashed line and filled circles). For £ = 0.5, 
even at S = 86 the difference between the simulations and 
the analytical result is still close to 10%. The deviations 
are quite pronounced for smaller values of the coupling 
parameter. The opposite is true for £ = —0.5. Here, the 
simulation results are close to the strong-coupling ana- 
lytical limit in the whole range of 2 values. For 5 < 4 we 
can discern, see inset in Fig. [JJ weak coupling behavior 
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that starts off with the PB-predicted value D* — 1 in the 
limit of S — > (Eq. solid line), that later follows the 
PB plus second-order corrections line (dot-dashed line) 
and then rapidly approaches the strong-coupling result 
D* ~ 0.92 (dashed line). As already noted the differ- 
ences between strong and weak coupling in this case are 
marginal in the rescaled representation. 

The dependence of D* on £ in Fig. [8] complements the 
above observations. The bound-state separation diverges 
for £ = both in simulations as well as in the analyti- 
cal limits. Here we reproduce the simulation results at 
( = 0.5 and —0.5, which again clearly show the conver- 
gence to the SC result for ( = 0.5 and the crossover from 
the PB result (solid line) to the SC result (dashed line) 
for C = —0.5 upon increasing the coupling parameter. 
The quantitative agreement between the simulations and 
the analytical results in the two limiting cases of PB and 
SC is excellent. Note here again that the second-order 
WC expansion (dotted line in the inset) begins to fail at 
small separations as it deviates from the limiting PB re- 
sults as well as from the simulation data. Since the PB 
pressure is zero on the PB phase boundary line, the va- 
lidity of the second order WC correction can be assessed 
analytically only by performing a two-loop calculation 
which goes beyond the scope of this paper. 



IX. CONCLUSIONS 

To summarize, we have derived theoretical forms for 
the interaction pressure as well as the counterionic den- 
sity profiles of asymmetrically charged planar surfaces 
with neutralizing counterions in between. Based on 
the field-theoretical methods we analyzed two different 
regimes of weak and strong coupling as defined by the 
electrostatic coupling parameter S. The crossover be- 
tween these two regimes is studied via Monte-Carlo sim- 
ulations. 

For small values of S, the system is described very 
well by the weak coupling (WC) theory, that in the low- 
est (zeroth) order coincides with the mean-field Poisson- 
Boltzmann (PB) result. The second order of WC cor- 
responds to a first-order loop expansion and represents 
the contribution from correlated quadratic fluctuations 
around the mean-field or saddle-point solution. This 
second-order correction, which is proportional to S, al- 
ways lowers the interaction pressure between the surfaces 
and thus leads to an attractive contribution to the total 
interaction pressure. Since it corresponds to an expan- 
sion of the partition function around the mean-field sad- 
dle point it has to be smaller, in absolute terms, than 
the mean-field result. Net attraction given by second- 
order fluctuations term (for interacting surfaces of equal 
sign) is therefore inconsistent with the nature of the WC 
approximation. 

For large values of the coupling parameter S, the weak 
coupling approach breaks down and the virial expansion 
amounting to the strong coupling (SC) approximation 



must be used. The SC approach is effectively a one- 
particle theory and takes properly into account the strong 
correlation and interaction of the counterions with ex- 
ternal surface fields on the leading order [9[. Following 
standard procedures, we derived an analytical expression 
for interaction pressure in the SC limit. The interaction 
pressure in this case is always lower than the PB result 
and can be negative (corresponding to a net attractive 
force) even for charged surfaces of equal sign. 

We compared both our theories, i.e. WC with second- 
order corrections and SC, with Monte-Carlo simulations. 
We found very good agreement for both theories in their 
expected regime of validity. As expected, the WC ap- 
proach is valid for separations D ^> 3. The second-order 
correction improves the small discrepancies between PB 
and MC results at large separations but it tends to fail 
for smaller distances when discrepancies get more pro- 
nounced. On the other hand, the SC theory describes 
the behavior perfectly at small separations D -C \/3. 
For small enough coupling parameter S, the validity of 
the PB approximation spreads to smaller separations, D, 
where PB and SC results nearly coincide (in the rescaled 
representation). Therefore, we may conclude that for 
sufficiently small S, the PB result is valid on the whole 
interval D. 

Note that the second-order WC correction term con- 
sistently diverges (toward large negative values) for small 
inter-surface separations, D < 1, irrespective of £, which 
makes it in general inapplicable in this limit. The rea- 
son for this is simple. For small inter-surface separations 
the mean-field solution becomes more and more homoge- 
neous, almost a constant, and the interaction free energy 
approaches its standard zero-frequency van der Waals 
form that diverges for small separations. 

In the case of surface charges with equal sign, £ > 0, 
the WC theory predicts no attraction and hence no 
bound state. The attraction and the corresponding 
bound state appear only for coupling parameters 5 that 
are large enough, as predicted by the SC theory. In the 
case of charged surfaces of opposite sign, ( < 0, the at- 
traction appears also in WC limit above a threshold value 
D* that represents the equilibrium surface-surface sepa- 
ration. 

It is notable that for charged surfaces of opposite sign, 
the WC analysis in general performs much better than 
for the surfaces of equal sign and that the SC and the 
WC results are very close to one another for charged 
surfaces of opposite sign in the whole range of rescaled 
separations, whereas they differ significantly in the case 
of surfaces of equal sign. There is also only a marginal 
difference between the PB and SC counterion density pro- 
files in the rescaled representation for surfaces of opposite 
sign. A reasonable explanation for this would be in our 
opinion that for oppositely charged surfaces the counte- 
rions mostly feel the effect of the strong uniform external 
field provided by the surface charges, which acts simi- 
larly in the strong as well as the weak coupling limit. 
Thus the mean-field and the strong-coupling approaches 
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should converge. In the case of similarly charged sur- 
faces, the mean-field theory depends more on the local 
counterion density whereas the strongly coupled counte- 
rions still feel mostly the external field. Thus the dif- 
ference between the WC and the SC frameworks in the 
C < and £ > cases. 

Our results support an emerging new paradigm, ac- 
cording to which the WC and the SC limit bracket the 
exact results for the interaction pressure between charged 
surfaces neutralized by mobile countcrions. They indeed 
provide quantitatively correct results for the interaction 
pressure in the limit of small and large inter-surface sep- 
arations, while at intermeditae separations the exact re- 
sults are always located between the two limits. It thus 
seems advisable that in analyzing the electrostatic inter- 
actions in colloidal systems one always calculates both 
analytic limits, the WC as well as the SC, in order to get 



a good handle on the range of values that the interaction 
can assume for any value of the electrostatic coupling pa- 
rameter. In future we intend to study the same system in 
the presence of added salt and dielectric discontinuities. 
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